% This function computes the obejctive function for estimation of Theta1 in
% the SR approach

function [Q,output,errorMes] = objectFunc(theta1InputValues,setup)
% errorMes = 0, everything is ok, 1 for reporting problems
if ~isfield('optimizer',setup)
    setup.optimizer = 2;
end

% Check if the lower and upper bounds are violated
if any(theta1InputValues < setup.lowerBoundsValues) || any(theta1InputValues > setup.upperBoundsValues)
    if setup.optimizer == 1
        Q = 1e6*NaN;
    else
        Q = 1D35;
    end
    output = [];
    errorMes  = 1;
    return;
end

% We turn theta1InputValues into a struct 
for i=1:size(setup.selectTheta1,1)
    theta1Input.(setup.selectTheta1{i,1}) = theta1InputValues(i,1);
end

% Accounting for calibrated coefficients, we unfold theta1Input by to get theta1 
theta1 = theta1Input;
% The setdiff implies that we do not overwrite elements in params based on calibrateparams
if ~isempty(setup.calibrateTheta1)
    namesCalibrate = setdiff(fieldnames(setup.calibrateTheta1),setup.selectTheta1);
    for i=1:size(namesCalibrate,1)
        name = namesCalibrate(i);
        theta1.(name{1}) = setup.calibrateTheta1.(name{1});
    end
end

% For the fully flexible model: Additional constraint on PHI to ensure identification
% The diagonal elements PHI must be increasing
errorMes = 0;
% for i=setup.nm+1:setup.nm*2+setup.nn-1
%     if theta1.(['phi',num2str(i),num2str(i)]) > theta1.(['phi',num2str(i+1),num2str(i+1)])
%         errorMes = max(errorMes,1);
%     end
%     % The Jordan form: normalization with (near identical eigenvalues under Q)
%     theta1.theta1.(['phi',num2str(i),num2str(i+1)]) = ...
%         (1-abs(theta1.(['phi',num2str(i),num2str(i)])-theta1.(['phi',num2str(i+1),num2str(i+1)])))^5000;    
% end
if errorMes == 1
    if setup.optimizer == 1
        Q = NaN;
    else
        Q = 1D35;
    end
    output = [];
    return;
end


% Solve the model
[alpha,beta,psi,mu,phi,sigmaxtil] = setQTSMparameters(theta1,setup);
if setup.modelNum == 1 
    %QTSM without interest rate smoothing
    [model,errorMes]  = solveQTSM(alpha,beta,psi,mu,phi,sigmaxtil,setup);
    %[model,errorMes] = solveQTSM_Realdon(alpha,beta,psi,mu,phi,sigmaxtil,setup);
    %[max(abs(modelR.g0-model.g0)),max(abs(modelR.gx(:)-model.gx(:))),max(abs(modelR.gxx(:)-model.gxx(:))) ]
elseif setup.modelNum == 2
    %QTSM with constant interest rate smoothing
    [model,errorMes] = solveQTSM_constRsmooth(alpha,theta1.rhor,beta,psi,mu,phi,sigmaxtil,setup);
end
if errorMes == 1
    if setup.optimizer == 1
        Q = 1e6*NaN;
    else
        Q = 1D35;
    end
    output = [];
    return;
end
model.modelNum = setup.modelNum;
model.macros   = setup.macros;
model.nm       = setup.nm;
model.nn       = setup.nn;
model.rLag     = setup.rLag;
model.factorSignResOn = setup.factorSignResOn;
model.rhor     = theta1.rhor;

% The Non-linear regression filer
T           = size(setup.yields,1);     % lenght of time series
ny          = size(setup.yields,2);     % max number of yields
nx2         = setup.nm + setup.nn;      % number of states for filtering
x2Hat       = zeros(nx2,T);             % array for storing estimated states
x2start     = ones(nx2,1)*0.01;         % the state at time t = 0 (only needed for NLS)
numIter     = zeros(T,1);               % number of iterations in the regressions
Qterms      = zeros(T,1);               % contribution to the objective function at period t
diffQValues = zeros(T,1);               % change in Q_t after the regression
residuals   = zeros(ny,T);              % the residuals
yHat        = zeros(ny,T);              % fitted bond yields
nyIndex     = zeros(T,1);
for t=1:T
    [x2Hat(:,t),Qterms(t,1),diffQValues(t,1),numIter(t,1),residuals(:,t),yHat(:,t),numObsPeriod] = ...
        NLSfilter(setup.yields(t,:)',model,x2start,setup.lambda,setup.tolFunNLSfactors,t);
    x2start = x2Hat(:,t);
    if model.factorSignResOn == 2
        % Penalty for getting negative weigths in the Taylor rule
        weightsTaylor = x2Hat(1:setup.nm,t);
        residuals(:,t) = residuals(:,t) + sum(weightsTaylor<0);
        Qterms(t,1) = residuals(:,t)'*residuals(:,t);           
    end
    nyIndex(t,1) = numObsPeriod;
end

% The average number of iterations in the regreesions
avgNumIter = mean(numIter,1);

% For output
Q                       = sqrt(sum(Qterms,1)/sum(nyIndex))*100;   %The objective function
yHatTrans               = yHat';
output.Qterms           = Qterms;
output.yHatTrans        = yHatTrans;
output.residuals        = residuals;
output.model            = model;
output.nyIndex          = nyIndex;
output.factors          = x2Hat;
output.macros           = setup.macros;
output.econAct          = setup.econAct;
output.avgNumIter       = avgNumIter;
output.numBootStep1And3 = setup.numBootStep1And3;
output.factorSignResOn  = setup.factorSignResOn;
%printToScreen(theta1InputValues);
end

